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ABSTRACT 

We present a statistical technique which can be used to detect the presence and prop- 
erties of moving sources contributing to a diffuse background. The method is a gen- 
eralization of the 2-point correlation function to include temporal as well as spatial 
information. We develop a formalism which allows for a derivation of the spacetime 
2-point function in terms of the properties of the contributing sources. We test this 
technique in simulated sky maps, and demonstrate its robustness in identifying the 
presence of moving and stationary sources. Applications of this formalism to the dif- 
fuse gamma-ray background include searches for solar system bodies, fast moving 
primordial black holes, and dense cores of dark matter proto-halos in the solar neigh- 
borhood. Other applications include detecting the contribution of energetic neutrinos 
originating in the solar system, as well as probing compact objects in long-timeline 
lensing experiments. 

Key words: methods: statistical — surveys — minor planets, asteroids: general — 
Kuiper belt: general — gamma rays: diffuse background — dark matter 



1 INTRODUCTION 

Diffuse background light is very important in understand- 
ing conditions and classes of objects in the Universe. This 
is due to the fact that the spectral, spatial, and amplitude 
information in a diffuse background is linked to the prop- 
erties of the otherwise unresolved contributing sources. For 
example, microwave background measurements include con- 
tributions of cosmic origin (Komatsu et al. 2009), as well as 
foregrounds of Galactic origin (Bernardi et al. 2005; Cruz 
et al. 2011; Dobler & Finkbeiner 2008; Dobler, Draine & 
Finkbeiner 2009; Gold et al. 2009). As another example, 
7-ray background measurements include contributions from 
unresolved blazars (Padovani et al. 1993; Stecker & Salamon 
1996; Miicke & Pohl 2000; Chiang & Mukhcrjcc 1998; Naru- 
moto & Totani 2006; Bhattacharya, Sreekumar & Mukherjee 
2009; Venters 2010; Ando & PavUdou 2009; Inoue & Totani 
2009; Dodelson et al. 2009), inverse Compton scattering of 
CMB photons by electrons accelerated at shocks around 
galaxy clusters and cosmic filaments (Loeb & Waxman 2000; 
Miniati 2002; Colafrancesco & Blasi 1998; Miniati, Koushi- 
appas & Di Matteo 2007), starburst galaxies (Pavlidou & 
Fields 2002), cosmic ray interactions with atomic and molec- 
ular gas in the Milky Way (Dermer 1986; Abdo et al. 2009), 
as well as the possible annihilation of dark matter (Ando 
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2009; Ando et al. 2007; Cuoco et al. 2007; Cuoco et al. 2008; 
Fornasa et al. 2009; Hooper & Serpico 2007; Lee, Ando & 
Kamionkowski 2009; Siegal-Gaskins 2010; Siegal-Gaskins & 
PavUdou 2009; Siegal-Gaskins 2008; Taoso et al. 2009; Bax- 
ter et al. 2010). 

Background events may be divided into two classes. 
Some events are generated by localized sources while oth- 
ers are generated by mechanisms which cannot be localized. 
In the first class the sources can be either spatially fixed 
(in celestial coordinates) or may exhibit proper motion (i.e. 
over a period of time their displacements are larger than the 
angular resolution of the detector). 

Using again the diffuse 7-ray background as an example, 
unresolved blazars, starburst galaxies, and emission from 
structure formation shocks would be considered spatially 
fixed sources of background. Cosmic ray events with inter- 
stellar gas would be considered a non-localized random pro- 
cess. Sources of background which will exhibit proper motion 
include the interaction of energetic cosmic rays with solar 
system bodies (e.g., small objects in the asteroid belt or ob- 
jects in the Kuiper belt and the Oort cloud) (Morris 1984; 
Moskalenko et al. 2008; Moskalenko & Porter 2009), dark 
matter annihilation around primordial black holes (Mack, 
Ostriker & Ricotti 2007; Lacki & Beacom 2010), and poten- 
tially nearby remnants of high density dark matter density 
peaks (Kousluappas 2006; Fieri, Bcrtone & Branchini 2008; 
Ando et al. 2008; Koushiappas 2009). In all these cases, in- 
dividual emission from any single object is not distinguish- 
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able, but the sum of these contributions may contribute to 
the diffuse 7-ray background. 

Correlations between individual events can help disen- 
tangle the contribution of various sources to the background. 
In this manuscript we present a formalism and a technique 
that can be used to identify the presence of background 
sources that exhibit spatial motion. In Sec. 2 we present an 
overview of the problem. In Sec. 3 we detail definitions that 
are used in the statistical techniques that follow. This allows 
us to write down the formal definition of the spacetime 2- 
point correlation function, which can be used to extract the 
moving signal in the diffuse background. In Sec. 4 we derive 
the form of the spacetime correlation function in 2 dimen- 
sions. A quantitive account of the uncertainty in the method 
is found in Sec. 5. In Sec. 7 we demonstrate the method's 
robustness in toy experiments and comment about the use 
of an instrumental point spread function. We generalize the 
formalism to realistic problems in 3 dimensions in Sec. 8, dis- 
cuss generalizations of the formalism in Sec. 9 and discuss 
applications and conclude in Sec. 10. 



2 OVERVIEW OF THE PROBLEM 

Suppose we have some objects moving on a 2-dimensional 
surface, each with a constant velocity. Every so often the ob- 
jects emit photons, which, when detected, we call "events". 
We record the location and time of each photon detection. 
The problem we are interested in is to take this collection of 
events and extract information about the objects: their exis- 
tence, their velocity distribution, their density distribution, 
and their event rate or luminosity (i.e. the rate at which 
photons are emitted from each object). 

There are two natural variables in this problem which 
ought to determine how difficult it will be to extract this 
information: the event rate of the objects and their number 
density. If there arc very few objects and their luminosities 
are very high it should be easy to identify the path of each 
object individually. In the opposite limit the objects' lumi- 
nosities are small but their number density is large. In this 
case it will be difficult to identify the sequences of events 
that trace the paths corresponding to individual objects. 

These two limits are represented in Fig. 1. Each panel 
in Fig. 1 is a plot of the location of all events in 10 arbitrary 
units of time^ . The left panel contains 5 objects each having 
a luminosity of 10 and an average speed of 5. The middle 
panel contains 50,000 objects each with an event rate of 0.01 
and drawn from the same velocity distribution as before. The 
right panel contains the same number of events as the middle 
panel, but they occur at random positions and times (i.e., 
there are no "objects"). In the left panel it is easy to measure 
the speed and event rate of every object (each generating 
about 100 events in total). This task is impossible, by eye, 
for the middle panel where each object generates 0.1 events 
on average. Indeed, it is even difficult to say whether or not 
the events come from objects at all, or if they are simply 

^ In these examples time and distance have arbitrary units and 
from now on these units will be set equal to 1. A phrase like 
"luminosity equal to 10" means an event rate of 10 per unit time; 
"an average speed of 5" means 5 units of distance per unit time, 
etc. 



generated randomly as in the right panel. In practice, the 
left panel is analogous to resolvable sources in the absence of 
any contaminating backgrounds while the middle and right 
panels represent diffuse backgrounds in the sky. Our goal 
is to be able to distinguish between the middle and right 
panels while learning something about the objects in the 
middle panel. 

The technique we employ is an application of the 2- 
point correlation function. One takes every pair of events 
and calculates their time separation and "velocity separa- 
tion" (their spatial separation divided by their temporal 
separation). One can then bin this data and make a 3- 
dimensional plot of number of pairs as a function of both 
time separation and velocity separation. The shape of this 
surface reveals information about the contributing objects. 
For instance, if all the objects are moving exactly at speed 
V, there will be lots of pairs of events whose velocity sepa- 
ration is V. The effect will be a ridge in this 2-dimcnsional 
parameter space. 

The situation can be made more realistic. Instead of 
the moving objects all having speed v, their speeds could 
be drawn from a distribution. Their event rates could also 
be drawn from a distribution. In fact wc might have many 
different populations of objects each having a different set 
of distributions for speed and luminosity. On top of this 
we could add a set of completely random events: a Poisson 
process such that there is some constant probability that 
an event occurs in any small region of spacetime. Below 
we will systematically discuss all these possibilities. First 
we present the simple 2-dimensional case with one class of 
moving objects along with a component of random events. 
This is the easiest way to present our formalism. Then we 
straightforwardly generalize to a realistic case where a dif- 
fuse background is made up of signals coming from various 
populations of objects as well as random processes. 



3 DEFINITIONS 

The analysis takes place on a 2-dimcnsional sky map, which 
is a collection of discrete signals that we define as "events" . 
Each event is assigned a spatial coordinate (position) and 
a time coordinate. For example, in the case of the Fermi 
Gamma-ray Space Telescope (Fermi), discrete signals are 
7-ray events recorded by the Large Area Telescope (LAT). 
The position is the location on the sky where the photon 
originated, and the time is the time of detection. It is im- 
portant to note that in realistic experiments the data comes 
not as a list of (position, time) for each event but as a list 
of (point spread function, time) for each event. The analysis 
that follows can be reworked for this more realistic situa- 
tion. However, we will start out by assuming that we simply 
have a collection of events where each event is specified by 
a position and a time. 

As wc arc interested in sources of events that can have 
velocities we also need a notion of distance. For realistic sky 
maps, the distance between two events is defined to be their 
angular separation. In our toy model with objects moving on 
a 2-dimensional surface, the distance between events is their 
Euclidean distance. We also define the "velocity separation" 
between two events to be the distance between them divided 
by their time separation. With these definitions, the appro- 
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Figure 1. Illustration of the two limits in the problem. The first figure contains 5 objects each with event rate 10 and the second contains 
50,000 objects with event rate 0.01. The third figure contains the same number of events as the second but they are distributed randomly. 
Naively, it is impossible to tell which of the last two figures contains random events and which contains moving objects. 



priate way to visualize the data is in a spacetime diagram 
where each event has both position and time coordinates. 

We will employ the 2-point function in a similar way 
to its use in galaxy-galaxy correlation studies. The galaxies 
correspond to what we have called events. To calculate the 
galaxy 2-point function for a particular angular separation 
9 one counts the number of pairs of galaxies in the sky map 
whose angular separation is between 6 and 9 ~\- A6'. That is, 
for every galaxy one looks in a ring of radius 9 and width A6 
around the galaxy and counts the number of other galaxies 
in this ring. The count is denoted by C{9,9 + A9){p), where 
p is an index labeling the central galaxy (see Fig. 2). If the 
events were distributed randomly one expects to find on 
average pV{9,9 + AS) galaxies in this ring, where p is the 
overall density of galaxies (number of galaxies in the sky 
map divided by the area of the map) and V{9, 9 + AS) is 
the area of the ring, equal to 27r(cos(6') — cos{9 + AS)). One 
then computes the correlation function ^ at separation 9 
according to 

.(n n , A .^ _ / C{9, 9 + AS)(p) - pV{9, S + AS) \ 

^(S,S + AS)-^ _____ ^, (1) 

where the average is taken over the index p of each galaxy. 
The correlation function ^{9,9 + AS) is interpreted as the 
fractional increase in probability (above random) that there 
is a galaxy in a ring between S and S -I- AS around any given 
galaxy. This is most easily seen by rearranging (1) into the 
form C = pV{l + ^). Notice that the correlation function 
is inherently a function of the shape and size of the ring in 
which the search for pairs of events is performed. 

Now we apply the 2-point function in our situation. We 
denote spacetime by 5 and we label spacetime events with 
the abstract index p, which carries all the information we 
have about the event. For example, for the event p, p{t) is 
the time the event occurred, p{x) is the x-coordinate of the 
event, etc. We define the spacetime 2-point function as fol- 
lows. For an event at p, let V{p) C S denote some volume of 
spacetime which is analogous to the shaded region in Fig. 2. 
When there is no confusion V{p) may also refer to the space- 
time volume of the region V{p). Two choices for V{p) are 
illustrated in Fig. 3. Let the number of events that occur 
within the region V{p) be denoted by C{p). When it is im- 




Figure 2. For a galaxy-galaxy correlation function we look in 
rings of a certain size centered on each galaxy and count the 
number of galaxies that lie inside each ring. The ring shown is 
V{p), centered on the galaxy (represented by the black x) having 
coordinates p. 

portant to remember that C{p) depends on the region V{p) 
we will write it as C{p; V). The spacetime 2-point function 
is then given by 

where the average^ is taken over every event in the sky map 

^ In order to be thorough we should really define ^ by ^(p; V) = 
{[C{p; V) — pV{p)]/ pV{p))i! , where the average is taken over an 
ensemble of Universes. Then we assume that our physical situ- 
ation is spacetime translation invariant so that ^(p; V) actually 
does not depend on the location p. Finally, in order to estimate 
^ from a set of data we claim that the average of ^(p; V) over an 
ensemble of Universes is equal to the average taken over all the 
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(i.e. over p) and V{p) denotes the spacetime volume of the 
region V{p). As before, p is equal to the overall spacetime 
density of events (total number of events divided by the 

spacetime volume of the sky map) . In a realistic application 
p will have dimensions of flux: events per square degree per 
time. 

If the events were all generated by a completely ran- 
dom Poisson process we would expect C{p; V) — pV{p) on 
average and (,{V) would bo 0. The 2-point function (,{V) 
is therefore to be interpreted as the fractional probability 
above random that the region V{p) contains an event given 
that there is an event at p. In the rest of this paper we will 
develop a formalism for deriving ^{V). 



4 2-DIMENSIONAL MODEL 
4.1 Ingredients 

Consider objects moving over a two dimensional surface with 
constant speeds and each having the same event rate (a 
"blinking rate" , so to speak) . Each event is then associated 
with an x, y, and t value and our "sky map" consists of the 
list of {x, y, t) for each event. The "blinking" of an object 
is a Poisson process with mean rate A: during the time dt 
each object has a Xdt chance of generating an event. Let the 
average density of objects be given by n, which has units of 
objects per area. The objects have speeds drawn from the 
distribution Pv{v): the probability for any given object to 
have a speed between v and v + dv is Pv{v)dv. We consider 
the case where the velocity distribution is isotropic (accom- 
modating the more general case, P^{v)d^v. is straightfor- 
ward). Finally, some fraction of the events will come from 
a random (Poisson) component with spacetime density po: 
there is a po dx dy dt probability of having such an event in 
any spacetime volume dx dy dt. 



4.2 The form of V{p) in 2 dimensions 

There are many possible choices for the spacetime region 
V{p). The simplest one is 



V{ri,r2; ti,t2)(p) 



{p eS: ti ^p'{t)-p{t) < t2 

A n d(p',p) < ra}, (3) 



where d(p ,p) is the spatial separation of spacetime events p 
and p' and A is the logical AND operator. This volume cor- 
responds to all the events whose temporal separation from p 
is between t\ and t2 and whose spatial separation is between 
ri and r2 a ring in spacetime with rectangular cross section 
(see left panel of Fig. 3). The volume of such a region is 
simply 

F(ri, ra; ti, fe) = 7r(r2 - r?)(t2 - ti). (4) 
A more convenient choice for V{p) is 
V{vuV2;ti,t2){p) = {p €S -.ti ^p'{t)-p{t) <t2 



events in our dataset. These are exactly the assumptions which 
must be made in the theory of galaxy n-point functions (referred 
to as ergodic conditions) . We will have more to say on the subject 
of estimators below. 



d(p'.p) 

\p'{t)-p(t)\ 



(5) 



This region is interpreted as the volume of spacetime that 
an object might explore between time ti and ta if it started 
at p and had any speed in the range from vi to V2- 
V{vi,V2; ti,t2){p) is illustrated in the right panel of Fig. 3 
for the case where p is at the origin of spacetime coordinates. 
If the event was an object which had velocity between vi and 
V2 its worldline would lie between the cones x^ + = vit^ 
and x'^ + y^ = V2t^ and we only consider the region where 
ti ^ \p'{t) — p{t)\ < t2. The volume of this region may be 
found by slicing the shaded region in the x — t plane, and 
rotating each small piece around the i— axis. The result is 



t2 V2t 



V{vi,V2; ti,t2) 



^11 



2-nxdxdt = -{vl - vi){tl -tl). (6) 



4.3 Coordinate systems 

The forms of V{p) presented in the previous section are eas- 
iest to visualize in the cartesian coordinate system [x, y, t) 
(see Fig 3). However, a more appropriate choice of spacetime 
coordinates are {v,t,(f)), defined by 



V = ^/{x^ + y')/\t\ 
t = t 

(j) = taia~^{y/x). 



(7) 



These arc just cylindrical coordinates (r, z, cf)) but with r 
scaled by the absolute value of z. The Jacobian for this 
change of variables is 

dV{v, t, (j)) = dx dy dt = vt^dv dt dcj). (8) 
The volume V of any region of spacetime V(p) is given by 



V= j dV{v,t,ct>)= j vfdvdtdct). 

V(p) V(v) 

For example, we can recover (4) as 

V{rx,r2\tut2)= / / / dV(v,t,<^) 

Jti Jri/t Jo 

and ( 6) as 

/'*2 /'^2 /'27r 

V{vi,V2;ti,t2)= I I I dV{v,t,4>). 



(9) 



Jti Jvi Jo 

For later use we define the {v, t, (f})p coordinate system 
which is the same as the coordinate system described in (7) 
except the center of coordinates {v = 0, t ^ 0) is at the 
spacetime point p. We also define the corresponding volume 
element dVp{v, t, (p) or dVp for short. The region dVp{v, t, (f>) 
is the infinitesimal version of (5), i.e dVp{v, t, 4>) = V{v, v + 
dv; t,t + dt), only not rotated about the t-axis. 

Finally we should note that it is just as easy to derive 
our results for rectangular coordinates. One just uses the 
coordinate system (vx,Vy,t) where 

Vx = x/t 

Vy = y/t 

t = t, (10) 
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Figure 3. The spacetime regions V{ri,T2\ ti,t'2){p) and V{vi,V2; ti,t2){p) where p is at the origin. The vertical axis is time and the 
horizontal axes are the x and y coordinates. In the left-hand figure the region between the two cylinders contains all events which have 
a radial distance from p between ri and r2. In the right-hand figure the region between the two cones represents the possible worldlines 
of an object starting at p and having a speed between vi and V2. Imposing a time separation between ti and t2 gives the filled regions. 



instead of (7). The analogue of (8) is then 

dV{vx, Vy,t) = dx dy dt — t^dvx dvy dt. (11) 



4.4 The form of S,{V) in 2 dimensions 

Our goal is to write down an analytic expression for (2). 
In the previous subsection we showed how to calculate the 
volume V{p). Now we move on to C{p;V), which can be 
interpreted as follows. Given an event at p, C'{p;V) is the 
probability^ of finding another event in the spacetime region 
V{p). 

There are two processes by which an event might oc- 
cur in V{p). Accordingly, we can break up C{p; V) into the 
sum of two terms: C{p; V) = [the probability of getting an 
event from an object that was at p] + [the probability of 
getting an event from any other source]. The first term can 
be thought about in a series of steps: given an event at p find 
the probability that it came from an object, that this object 
moves into the region V{p), and that this object triggers a 
new event while in this region. 

The probability pi that any given event came from a 
moving object (as opposed to being generated by the Poisson 
component of the background) is the ratio of the flux from 
moving objects to the total flux: 



nX 



Pi 



Pi 
P '' 



nX + po 

where pi = nX is the average flux of the moving objects. As 
before, p is the total spacetime density of all events (i.e. the 
overall flux). 

The probability that the object moves into the region 
dVp{v, t, (j)) is simply the probability that its speed is between 
V and V + dv, Pv{v)dv, multiplied by d0/27r, the probability 



^ Or, if it is greater than 1, C{p;V) is the expected number of 
events in the region V{p). 



that it is moving in a direction between cf> and (j) + d(fi ■ li 
the object makes it into the region V{p) the probability of it 
generating a second event is Xdt. Therefore, the probability 
that there was an object at p which moved into V{p) and 
generated another event is^. 



Pi 



Pv{v)dv Xdt 



d4 
2^' 



(13) 



V(p) 



The second term in C{p; V) is simply pV{p), where p dx dy dt 
is the probability that any random spacetime volume 
dxdydt contains an event from either an object or the ran- 
dom component (note that p — pi + po). 

Therefore, putting together these parts and plugging 
them into (2) we find 



m = 



/ C{p; V)-pV{p) \ 
\ PV{P) I 
Pi P.,{v)dvXdtd^l2-K 



pV{p) 
P^{v)dvXdt d(p/2Tv 



(pi + Por Vip) 



(14) 



As is usually done for galaxy-galaxy correlation func- 
tions let's see what happens when we take the limit V{p) 
dVp{v,t,(j}). Using (8) we see that 



aV)^adV{v,t,c^)] = 



Pi A P^v) 



2n{pi + Po) 



yt2 



(15) 



^ If the objects do not have an isotropic velocity distribution then 
this probability is Pi;{v, (p)dvd(l>, where Pj(w, </>) is the probability 
density for the velocity vector. 
^ In full generality this 

where 



equation would be 
X(t)dt is the probabil- 
ity that an object which generated an event at t = generates 
another event in the time interval between t and t + dt. 



Pi /^(pj V^(v,4,)\(t)dvdtd(j>, 
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This limit is finite and the function ^ traces the veloc- 
ity distribution of the population of objects. Therefore, if ^ 
can be measured for multiple values of v then it is possible 

to directly reconstruct both the velocity distribution of the 
moving objects and information about their abundance and 
luminosity. 



5 THE ERROR IN ^ 

Getting a handle on the error in a measurement of ^(V) 

is just as important as calculating ^{V) itself: any practical 
application of this method will reveal nothing if the uncer- 
tainty in £,{V) is comparable to ^{V). The zeroth order dis- 
covery that can be made using the 2-point function is the 
detection of the presence of moving objects. This is done 
by rejecting the hypothesis that (,{V) = for all choices of 
V{p), which is possible oidy if A^/^(V) < 1 for some choices 
of V{p). An estimate of the error is also essential when fitting 
the theoretical value for ^ to the data; i.e. when performing 
a fit to determine the physical parameters describing the 
density, luminosity, and velocity distribution of contributing 
sources. 

Fortunately, the errors in correlation functions have 
been thoroughly studied in the case of galaxy-galaxy cor- 
relations (Peebles 1973; Fall & Tremaine 1977; Mo, Jing 

6 Boerner 1992; Fry & Gaztariaga 1993; Landy & Szalay 
1993; Hamilton 1993; Bernstein 1994; Szapudi & Colombi 
1996; Szapudi & Szalay 1998). We emphasize that all the 
technology that has been developed for calculating 2-point 
functions for galaxies and quantifying their errors can (and 
should) be straightforwardly applied to our 2-point function. 
As stated before, the only conceptual difference between the 
two tools is the choice of V{p). 

In particular, we apply the results of Landy & Szalay 
(1993) (hereafter LS93) to the present problem. In the ex- 
amples below we measure (,{V) using an unbiased estimator 
which is identical to the DD/RR ratio in LS93. This is done 
for simplicity. The unbiased LS93 {DD - 2DR + RR)/RR 
estimator was shown to have a smaller variance and should 
be used in practical applications. (In LS93, DR refers to the 
cross-correlation of the observed data with sets of completely 
random data, while DD and RR axe the auto-correlation 
functions computed for the data and for a completely ran- 
dom set of data, respectively.) 

To quantify the error in i,{V) we adapt the LS93 expres- 
sion for the variance of the {DD — 2DR+RR)/RR estimator 
for small correlations (i.e. small values of (,{V), likely in cases 
of physical interest). For a given shape V{p) (with spacetime 
volume V) the variance of the estimator is given by 



NpV 



(16) 



where A'' is the total rmnibcr of events in the sky map. This 
can be seen to be the same as Eqs. 43 and 48 in LS93 by 
writing p = N/V, where V is the total spacetime volume 
of the sky map and noting that V{p)/V is equal to LS93's 
Gp{9). The signal to noise ratio is then 



(17) 



These expressions should be used to determine the optimal 



volumes V{p) for any given application. Ideally, V{p) should 
be chosen to make the signal to noise ratio large while keep- 
ing V{p) small enough that many choices for V{p) can be 
measured for the sky map. This dilemma occurs with galaxy- 
galaxy correlation studies as well. An annulus of specific size 
(sec Fig. 2) corresponds to one choice of V(p). One would 
like to choose the width of the annulus as small as possible 
so that the correlation function can be measured at many 
different angular scales. However, the smaller the width of 
the annulus the larger the uncertainty in the measured value 
of the 2-point function. 

As is well-known in galaxy-galaxy correlation studies, 
measurements of ^ at two different angular sizes can be 
highly correlated. This issue will also affect any measure- 
ment of the spacetime 2-point function: the measured ^{V) 
for different choices of V{p) may be correlated. Therefore, a 
fitting to extract physical parameters should include an 
estimate of the covariance of £,{V) between different V(p)'s. 
A variety of methods have been developed to estimate or 
predict this covariance matrix. Many of these are trivially 
adapted for use in this case. Bootstrapping (e.g. Ling, Bar- 
row & Frenk (1986); Barrow, Bhavsar & Sonoda (1984); 
Fisher et al. (1994)) and jackknife resampling (Lupton 1993; 
Zchavi ct al. 2002) require measuring the correlation func- 
tion on various subsets of the full data set and analyzing the 
variation among these estimates of ^. If generating fake data 
sets is feasible then one can simply measure the correlation 
function on many fake maps to find the covariance of ^(V) 
between various F(p)'s. 



6 POINT SPREAD FUNCTION AND 

COMPUTATIONAL CONSIDERATIONS 

In this section we discuss two ways to include information 
about the point spread function (PSF) into the derivation of 
the form of the 2-point function. This will serve as a guide 
for incorporating the PSF in realistic applications. 

The PSF of a detector quantifies the uncertainty in 
its measurement of the locations of events in spacetime 
(Morales, WilUams & De Young 2004). The PSF typically 
takes the form PSF(pt — Po), where pt is the true loca- 
tion of the event and po is the location that the detector 
reports, the "observed" location. The PSF is a probabil- 
ity density on spacetime: PSF(pt — Po)dV is the probability 
that if the detector reports an event at po it actually arrived 
from the spacetime region dV centered at pt^. As a result, 
J PSF(pt —Po)dVt = 1. Additionally, there is the probability 
Pd that if a signal (e.g. a photon) arrives at the detector it 
will actually be detected as an event. 

If we are given the PSF for a given event we can do a 
more precise job of computing C{p; V). As above we want 
to answer the question: given that the detector reported an 
event at po what is the probability that the detector reports 
another event in the spacetime region V{p)7 

If the detector reports an event at p there is a pi chance 
that it received a signal from a moving object. But the true 

^ Because the time resolution of detectors is generally excellent 
compared with the spatial (or angular) resolution, the PSF is 
usually given as a function of spatial coordinates only. The PSF 
we have defined would then be equal to <5(tt — to) PSF(rt — ^o). 
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location of the object could be anywhere, with probability 
given by the PSF. The object can have any velocity and can 
emit a signal at any later time. This signal has a pa chance of 

being detected. The location of the observed event is again 
determined by the PSF. Specifically, we have 



C{p;V) 



pi J PSF{pt-p)dVp{vt,tt,ct>t) 

Ptes 



)Xpddv' dt' d(j)' 



PSF{p't-po)dVj,,{vo,to, 



(18) 



Paev(p) 



In words, there is a pi chance that the observed event at p 
came from a moving object. Given that it came from a mov- 
ing object there is a PSF(pt — p) dVp{vt,tt,4>t) chance the 
event actually occurred in the region dVp{vt,tt,(l)t) around 
the point pt = {vt,tt, (f>t)p (recall the definition of dVp 
at the end of the section on the choice of V{p)). Then 
there is a {l/2'K)Pv{v')Xpddv' dt' d(j)' chance that the ob- 
ject moves into the region dVp^{v' ,t' ,4>') around the point 
p't = {v' ,t' ^(j))p^ and emits a signal which is reported by 
the detector. Finally there is a PSF(pJ — po) dVp' {vo,to, (j)o) 
chance that this event is reported as having occured in the 
region dVpi^{vo,to,4>o) around the point po = (wo, io, ^o)pj- 

All the possibilities are taken into account by integrat- 
ing pt and p[ over all of spacetime (the object could actu- 
ally have been located at any point and could have moved 
to any other point) and by integrating po over the region 
V{p) (we are only interested in the possibilities where the 
detector reports the second event in the region V{p)). For 
clarity we have omitted the pV{p) term in C{p;V), which 
represents the probability of a reported event in V{p) from 
any source besides an object moving from p into V{p). One 
can show that (18) reduces to the numerator of (14) when 
PSF(pt - Po) = S{pt - Po) and pd = 1. 

The spacetime correlation function is an example of a 2- 
point correlation function and so any method that is used to 
compute 2-point functions may also be used here. In galaxy- 
galaxy studies, the galaxies are localized sources and the 
2-point function is measured by counting pairs of galaxies 
which have a particular separation. When looking for mov- 
ing objects using gamma-ray data, for instance, the events 
are also localized. Computational procedures then carry over 
directly. Typically, counting pairs of events is an iV^ pro- 
ccss'^. For example, in gamma-ray diffuse studies the num- 
ber of events is proportional to the observation time as well 
as to the effective area of the detector. 

In other situations the data do not come as localized 
events but as a continuous amplitude across the sky. This 
case can be treated by first discretizing the survey area 
into small "cells" or pixels. Each pixel now has a contin- 
uous value. For a particular V{p), the correlation function 
is found by multiplying the value of the pixel p by the sum 
of the values in the pixels in the volume V(ji). The expected 



We point out that efficient algorithms with better than N'^ 
scaling have recently been developed. See, for example, Moore 
et al. (2001); Zhang & Pen (2005); Eriksen et al. (2004). 



value of this quantity (i.e. the denominator in (2)) is the 
average pixel amplitude squared multiplied by the volume 

of y(p). 

This method of computing the 2-point function can be 
used as an alternative way to account for the detector point 
spread function. Following Morales, Williams & De Young 
(2004), every discrete photon event in spacetime is replaced 
by its point spread function. The overlap of the point spread 
functions for all observed events forms a continuous density 
over the survey area and observation window. The 2-point 
correlation function for any choice of V{p) can then be mea- 
sured as described above. We note that this method suffers 
no performance penalty for increased numbers of observed 
events because the events are essentially binned into pixels 
in spacetime, with each pixel having a value given by the 
linear superposition of all contributing PSFs. 



7 EXAMPLES OF THE 2-DIMENSIONAL 
FORMALISM 

In this section we will demonstrate the accuracy of the 
derivations by measuring ^ for three different simulations in 
which the objects move according to a specific speed distri- 
bution. A generic choice for Pv is the Rayleigh distribution: 
the speed distribution for a 2-dimensional isotropic Gaussian 
velocity distribution. It has the form 



(19) 



with mean speed v = We choose V{p) to be the re- 

gion described by (5) and shown in the right panel of Fig. 3. 
We note that any choice of the shape of V{p) is allowed. 
The shape V{p) used here is adapted to the search for ob- 
jects which move in straight lines at constant speed. For 
sources with different patterns of motion, other choices for 
V{p) may be more appropriate. However, because the choice 
affects the counting of pairs of events when measuring ^ it 
nmst be taken into account in the theoretical derivation of 

With these choices the integral in (14) becomes 



Py{v)dvdt 



V(p) 



27r 



a-' 



dvdt 



27r 



tl VI 

= (t2-ti) e" 



Inserting this expression into (14) and using (6) we find 



^(i'i,i'2; ti,t2) = 



PlA {t2 - tl) 



v\l2ar 



2 /2a 



(7r/3)(pi+po)2(i;§-t;?)(ti-i?) 



(20) 



Given an event map we can measure C(«i, V2\ t\,t2) for 
any choice of the 4 parameters («i, W2, ii, ^2). In practice, 
a fit can be attempted in order to discover the 4 physical 
parameters A, pi, po, and a. While p = pi -|- po is measured 
directly the parameters pi and A are combined as a single 
normalization factor and so the most a fitting analysis would 
reveal would be the combination piA. In the 2-dimensional 
case this is true for any choice of V(j))-, as can be seen from 
(14). Of course, knowledge of any one of A,pi, or po can be 
used to find the other two. 



© 0000 RAS, MNRAS 000, 000-000 



8 Alex Geringer-Sameth and Savvas M. Koushiappas 
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Figure 4. A toy example demonstrating the use of the spacetimc 
correlation function to discover the presence of localized event 
sources with non-zero speeds. The t = slice of ^{v,t) is plot- 
ted showing the theoretical prediction (red x's), the measured 
value (blue squares), and the measured value for the case of com- 
pletely random events (black triangles). The hypothesis that the 
pattern of events in the sky map is Poisson {^{v, t) = 0) is clearly 
rejected at high significance. The error bars in the measured quan- 
tities are explained in the discussion surrounding (16). The sky 
map contained 3.5 million events, all from moving objects, though 
each object contributed only 0.1 events on average. The blue data 
points are measured from a larger version of the map shown in 
the central panel of Fig. 1 while the black points are measured 
from a larger version of the map shown in the right panel. 



mean speed w = 5. The objects then have the same density, 
event rate, and speeds as in the middle image of Fig. 1. 
The expected number of events triggered by each object is 
0.1 which means that although there are about 35 million 
objects present, less than 10% of them will trigger even a 
single event. Overall, there are roughly 3.5 million events in 
our sky map. 

In the measurement of ^{vi ,V2; ti,t2) we take 112=111-1- 
1 and t2 = ti + 1, i.e. we choose non-overlapping bins of size 
1 in both time and velocity separation. The subscripts on vi 
and ti are dropped and ^(wi,i'2; t\,t2) is relabeled ^{v,t). 
The 2-point function is then measured for « = 0, . . . , 19 and 
for t — 0,1, 2. The t — slice of the measured ^(«, t) is shown 
in blue in Fig. 4 along with the theoretical value (20), shown 
in red. The black curve is the 2-point function measured 
for a sky map containing the same number of events but 
placed randomly. The separation V2 — vi — 1 was selected 
for illustrative purposes. The time separation t2 — ti = 1 
was then chosen to be close to the optimal separation found 
by maximizing the signal to noise ratio (17) for ti = 0. 
The error bars are computed according to (16). This is a 
slight abuse since the estimator plotted is DD/RR and not 
{DD — 2DR + RR) / RR. In practice it is recommended to 
use the latter estimator. 

It is clear that moving objects are detected at a very 
high significance (i.e. the hypothesis (,{v,t) = is rejected). 
The measured value ^(0,0) — 0.15, for example, is about 15 
standard deviations from ^ = 0. A fit to recover the param- 
eters A, pi, and a can be attempted using ^{v,t), which is 
measured at the lattice of points {{v, t) : v = 0,1, . . . ;t = 
0, 1, . . .}. In practice, the full covariance matrix of errors be- 
tween different u-bins should be included in such a fit (see 
last paragraph in Sec. 5). 



Our first simulation will only contain moving objects. 
In the second we will add a component of random noise and 
in the third simulation both random noise and a popula- 
tion of stationary objects will be considered in addition to 
the moving objects. One of the goals of these simulations 
is to demonstrate that the 2-point function can tell the dif- 
ference between a background containing a population of 
sources and a completely random background. We do this 
by generating a second sky map for each example with the 
same number of events but distributed completely randomly 
throughout the spacetime volume. The 2-point function is 
measured for this randomly generated sky map and is plot- 
ted along with the 2-point function measured from the ac- 
tual simulation. If the events are randomly generated there 
should be no correlations at all: £,{vi,V2; t\,t2) should be 
for all values of 11, 12, fi, and t2. 

7.1 Example 1: Moving sources only 

In the first simulation there is no random noise: po = 0. 
We simulate an area with dimensions 13,200 x 13,200 for 
time 10. The density of objects is n = 0.2 and each has 
an event rate A = 0.01 yielding an estimated flux of pi — 
nX — 0.002 events per unit area per unit time. Their speeds 
are distributed according to a Rayleigh distribution with a 



7.2 Example 2: Moving sources and a random 
component 

Let us see if the spacetime 2-point function can tell the dif- 
ference between a collection of moving objects plus random 
noise and a situation with just random noise, where both 
cases have the same total flux p. The sky map has the same 
dimensions as before and the moving objects have the same 
number density, luminosity, and speed distribution as before 
yielding pi = 0.002. We choose the random component to 
have the same flux po = pi so that p = po + pi = 0.004. 
There are about 7 million events in this sky map, half com- 
ing from objects and the other half coming from the random 
component. 

In the calculation of (,{vi,V2; ti,t2) we choose V2 = vi-\- 
1 and t2 = ti -1-0.9. The results are plotted in the left panel of 
Fig. 5. Again it is clear that the moving objects are detected 
even in the presence of random signal in the sky map. 

How impressive is this result? Could we have just looked 
at the data by eye and spotted the presence of moving ob- 
jects? If each object generates at most a single event then 
clearly it is impossible to determine anything about their 
motion or to distinguish this from the case of completely 
random events. The fraction of events which come from ob- 
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Figure 5. The 2-point function ^ measured for two simulations color-coded as in Fig. 4. Each contained 7 million events. Objects had 
the same event rate as the first simulation. Left: Moving sources and random noise. Half the events came from moving objects and half 
were generated completely randomly to represent noise. Right: Moving sources, stationary sources and random noise. A third of the 
events are from moving objects, a third from stationary objects, and the last third were generated randomly. 



jects that trigger more than one event is 



P>i = 



nAY^k n{k;\T) 

k = 2 

(pi + pq)AT 

nA Ar(l-e~^^) 
(pi + po)AT 



pi + Po 



(1 



(21) 



where A total area of the sky map, T is the observation 
time, and 7r(i;M) = e~^' NP /i\ is the Poisson distribution 
with mean M. In our case, pi — po — 0.002, A — 0.01, and 
T = 10. Substituting these values into (21) gives P>i — 
0.048. That is, less than 5% of the events in our simulated 
sky map come from objects which generate more than one 
event. Furthermore, 95% of these events come from objects 
which generate exactly two events during the time T. If one 
was to try to spot individual moving objects in the sky map 
one would need to be able to take 200 events and out of the 
nearly 20,000 possible pairs of these events spot the 5 pairs 
which correspond to an object triggering an event, moving, 
and triggering a second event. 

We can illustrate this difficulty by examining a smaU 
area of the sky map from our simulation. The left panel of 
Fig. 6 shows all the events that occurred in a 150 x 150 area 
during the entire 10 units of time. On the right the same 
events are shown but are identified as having come from 
objects (blue) or as random events (orange). Events which 
came from the same object are connected with a line. The 
difficulty of discovering moving objects by eye is evident. 
The 2-point function is statistically able to pick up on the 
rare occurrences where an object generates more than one 
event. 



7.3 Example 3: Moving sources, fixed sources, 
and a random component 

As a final example, a third class of objects are added to 
the simulation. These are stationary objects which do not 
move during the course of the observation. The presence of 
such objects should manifest itself as a spike in the 2-point 
function at t; = 0. 

The dimensions of the sky map are the same as in the 
two previous examples. The moving and stationary objects 
have the same event rate, A — 0.01, and spatial density, 
n ~ 0.133. The moving objects have the same velocity dis- 
tribution as before. The random component has spacetime 
density po ~ 0.00133. Therefore, the total density of events 
is p = po -f pi + p2 = 0.004, which is the same as in the last 
example. The subscript 2 denotes the stationary objects. 
Each component contributes roughly the same number of 
events to the sky map. 

Including the stationary objects into the 2-point func- 
tion just requires replacing the Rayleigh distribution with 
the Dirac delta function centered at t; = 0; Pv{v) — S{v). 
Eq. 20 becomes. 



^2{V1,V2; ti,t2) = 



3p2A {t2 — tl)(?ui,0 



TT p 



(22) 



and (,i{vi,V2; ti,t2) is given by (20) except that the total 
density in the denominator includes the stationary objects, 
p — po + pi + P2- The function (5„i,o is if wi > and is 
1 if Hi = 0. The measured 2-point function ^(^1,^2; ^1,^2) 
is simply the sum of the 2-point functions for each class of 
objects: ^ = ^1 + ^2. 

As before we choose V2 — vi + 1 and t2 = ti + 1. The 
results are plotted in the right panel of Fig. 5. The spike at 
V = due to the stationary objects is apparent. Its height is 
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determined by both and ^2- Since the shape of ^ when v > 
can be measured the contribution of the moving objects 
to the spike aX v = can be subtracted. 



8 OBJECTS IN 3 DIMENSIONS 

The 2-dimensional situations examined so far are, of course, 
only toy models for astrophysical applications. In this sec- 
tion we develop a more realistic theory of the use of the 
2-point function. The derivation of the form of ^ is based 
on precisely the same arguments as in the 2-dimensional 
case. Simulations analogous to those in the previous section 
can also be performed in three dimensions and will agree 
with the theoretical form of ^. In performing an actual mea- 
surement of ^ simulations should be tailored to the specific 
application. We defer such detailed modeling to future work 
wherein we will apply the formalism to all-sky gamma-ray 
data Geringer-Sameth & Koushiappas (2011). 

A diffuse emission all-sky map (e.g. from Fermi-LAT) 
is a 3-dimensional representation of a 4-dimensional process 
since we cannot measure line-of-sight distances for individ- 
ual events. The distance to a source determines both its flux 
on Earth and its angular speed across the sky. This coupling 
between distance, speed, and flux is what makes the analy- 
sis more complicated. Previously, the velocity distribution 
Pv{v) and the luminosity A were independent quantities. 
Now we must consider probability distributions which de- 
pend on both II and A: closer objects have higher angular 
speeds and look brighter than distant objects. While the 
spacetime 2-point function in this situation is still defined 
by (2) it is more difficult to derive the analogue of (14). The 
analysis of this section will develop the theory of the space- 
time 2-point function in the case of a realistic sky survey. 



8.1 Summciry of the mesisurement of ^ 

The computation of ^ proceeds exactly as in the 2- 
dimensional case. The sky map consists of events, each hav- 
ing a directional coordinate (the apparent direction of the 

photon's origin) and a time coordinate. The "distance" be- 
tween events is the angle between them measured along a 
great circle. The velocity of interest is now an angular ve- 
locity: the "velocity separation" of two events is defined as 
the angle between the two events divided by their time sep- 
aration. The sky map is again a spacetime diagram, though 
not with the usual rectangular coordinates for the spatial 
axes. It can be visualized as a series of concentric spheres, 
each representing the celestial sphere, with different spheres 
corresponding to different slices of time (with t increasing 
as the radius of the spheres increases). In this picture the 
worldlines of objects moving at constant angular speed are 
Archimedean spirals in spacetime. The volume of a region 
in this spacetime has units of solid angle x time. 

In the next sections we derive expressions for V{p) and 
C (p; V) , the latter in terms of parameters describing the 
various populations of objects which contribute to the sky 
map. 



8.2 The form of V{p) in 3 dimensions 

One defines V{p) as some volume of spacetime <S. Convenient 
choices include 

V{eu02;ti,t2){p) = {p € S : h p' (t) - p{t) < t2 

A6'i ^ d(p',p) < 6I2}, (23) 

where p{t) is the time coordinate of the event p and d{p' , p) 
is the angle between spacetime points p and p', and 

V{ui,uJ2;ti,t2)ip) = {p £ S : ti p (t) - p{t) < t2 

where wi and 0J2 are angular speeds. These are the analogues 
of Eqs. 3 & 5. In (23), V{p) contains the events which occur 
in an ammlus around p with inner and outer radii 61 and 62 
and which occur in the time interval p{t) + ti to p{t) + t2- 
Note that when ii = and ^2 = 00 this region is exactly 
that used for galaxy-galaxy correlation studies. In (24), V{p) 
represents all the events which could have been triggered by 
an object moving from p if it had an angular speed between 
wi and W2 and with the same time separation constraint. 

If we are looking at a small area of the celestial sphere 
that can be approximated as fiat space then we can choose 
V(p) to be "anisotropic", i.e. choose volumes such as 

V{v:^i,Vx2\ Vyi,Vy2\ ti,t2) = {p £ S : (25) 

tl ^p'{t)-p{t) <t2, 

p'(x) — p(x) 

p'(t)~p(t) 

, „ ^ p'{y)-p{y) , 1 
^''y'^ p'^t)-p{t) 

This choice of Vijp) is useful when a class of moving objects 
has an anisotropic velocity distribution, or when the proper 
motion of the earth or the detector is important. 

The volume of the region V(p) is calculated in a way 
similar to the 2-dimensional case. For instance, the volume 
of the region specified by (23) is found by first computing the 
solid angle of the annulus between 9i and 82 , and multiplying 
this by the time interval: 

¥{01,92; UM) = 27r(cos6ii - COS6I2) (t2 - ti). (26) 
The volume specified in (24) is slightly more complicated: 

y(wi,W2;ti,t2) = / dt d(j> sinOde 

J ti Jo J ujit 

„ sin(wii2) — sin(a;iti) 

= ZTT 

27r sin(a;2t2) - sin(a;2ti) ^27) 

U!2 

Equations 26 and 27 hold only when 62 < and W2i2 < tt, 
respectively. Otherwise the annulus begins to overlap itself. 
This is only an issue if one is searching for objects which 
moved across the entire sky during the observation period. 

In the limit where t2 ^ t\ + dt and W2 — > wi -|- (27) 
becomes (dropping subscripts) 

y(aji,cj2; ti,t2) -)■ dT/(w,t) = 2-Ktsm{ojt)dujdt. (28) 
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Figure 6. Left: All events whieh oecurcd in an area of the sky map with dimensions 150 X 150 during the entire observation time. Right: 
The same events but identified as objects (blue) and random events (orange). Events which came from the same object are connected 
with a line. Less than 5% of events come from objects which caused more than one event. 



This is the analogue of the 2-dimensional (8). 

As in the 2-dimensional case we now define a convenient 
coordinate system for every point on the celestial sphere. 
The coordinates (oj, (j>, t)p are related to the global celestial 
coordinates (plus time) as follows. First we consider a ro- 
tated set of spherical coordinates (G, $)p in which p is at 
the north pole and the line $ = intersects the north celes- 
tial pole. That is, the new and old coordinates are related 
by a rotation in which p slides along a line of longitude to 
the north celestial pole. Then new coordinates {u),(f),t)p are 
related to (Q, $) by = u}t and ^ = <j>. This is a mapping 
from {uj,(j>,t)p to the global celestial coordinates (the time 
coordinate is unchanged) . Using (28) we can write down the 
volume element in these coordinates. The spacetime volume 
(solid angle x time) between u} and uj + du, between and 
(f) + d(j>, and between t and t + dt is 

dVp{(jj, cj>,t) = t smiujt) diu d(j> dt, (29) 
and one can check that (27) is recovered as the integral 

dVp{uj,cj>,t). 

8.3 Ingredients needed to derive C(p\V) in 3-D 

Besides the volume V{p) we need to derive an expression for 
C(p;V) in (2). This quantity depends on the properties of 
the sources which contribute events to the sky map. Gen- 
erally, the sky is populated by different classes of objects, 
each with its own velocity distribution, luminosity function, 
and spatial distribution. Let's denote the different classes of 
objects by the subscript i. Then for each class we define the 
following functions. 

• Pi^L{L) dL is the probability that an object of class i 



has an intrinsic luminosity between L and L + dL. L is the 
number of photons per second emitted by the object. The 
distribution Pi.L{L) is normalized to 1: Pi^L{L)dL — 1. 
The function P^^l is commonly called the luminosity func- 
tion of the population. 

• ni{R,fl) is the physical number density of i-type ob- 
jects which lie a distance R away from the detector in the 
direction Cl on the celestial sphere. This quantity has units 
[length]"^. 

• fi{v; Cl) specifies the tangential velocity distribution of 
i-type objects. The quantity fi{v; Cl)d^v is the probability 
than an object of class i located in the direction Cl on the ce- 
lestial sphere has its tangential velocity vector in the range 
d^v around v. These velocities are proper velocities, mea- 
sured relative to the earth (or the detector). "Tangential" 
means that the velocity is perpendicular to the line of sight*. 
This distribution is also normalized to 1. 

• p, defined above, is the average number of events de- 
tected per solid angle per time. It can be estimated from 
the sky map by dividing the total number of events by the 
time over which the sky map was measured and by the to- 
tal solid angle of the map. For maps with large numbers of 
counts this estimator will be adequate. In practice, one may 
need to modify the procedure for surveys with unequal ex- 
posures across the sky: the quantity p may be position and 
time-dependent. If we divide p by the detector area A we 
get p, the total flux per solid angle. 



" We have implicitly assumed that the line of sight velocity of 
any object is small enough that the change in its distance does 
not affect its flux. That is, the objects are all far enough away 
so that viosi/R <^ 1, where f is a measure of the time separation 
between the region V{p) and the event p. 
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8.4 Derivation of C(p; V) in 3 dimensions 

The quantity C(p; V) is tiie probability of finding an event 
in the region V{p) given that the detector reported the event 
p. It can also be thought of as the expected number of events 
in V{p), given an event p. First wc break C'{p;V) into the 
sum of 2 terms; C{p; V) = [the probability that the event p 
was caused by an object which moved into the region V{p) 
and triggered another event] + [the probability of finding an 
event in V{p) for any other reason]. As in the 2-dimensional 
case the second term is simply pV{p). 

The first term can be broken up into the product of 3 
probabilities: [Ci: the probability that the event p came from 
an object of class i with luminosity L located a distance R 
from the detector] x [C2: the probability that this object 
has a velocity that takes it into the region V{p)] x [C3: the 
probability it triggers an event while in V(p)]. The product 
then needs to be integrated over R and L and summed over 
i. 

The first factor, Ci, is the ratio of photons received 
from i-type objects with luminosity L and distance R in the 
direction tip to the total number of photons received from 
the same direction: 



Ci = ^ \ni{R,Clj,)R' dR] [Pi,L{L)dL] 



LA 



(30) 



L47r7?2_ 

where A is the effective area of the detector. It is worth 
noting that Ci does not actually depend on A since p will 
also be proportional to A. 

The second factor, C2 , is the probability that an i-type 
object will have a velocity which takes it into the region 
V{p). We will have to integrate over a range of velocities 
which correspond to the object moving into V(p). It will, 
therefore, be useful to use the coordinate system defined in 
the discussion leading to (29). We can adapt the velocity 
distribution ,fi(v; Cl) to the now coordinates by introducing 
the function fi{v, (f>; Cl) defined so that 



fi{v,4>; Cl) dv d(f> = fi{v; Cl)d v. 



(31) 



The quantity fi{v,<j); d)dvd4> is to be interpreted as the 
probability that an object of type i has tangential speed be- 
tween V and V + dv and is moving in a direction between 
(j) and (j) + d4>, where refers to the coordinate label in our 
new coordinate system whose north pole coincides with the 
direction Cl as described previously. Next we relate the dis- 
tance to the object to its angular velocity using the relation 
V = Rw. Thcroforo, the quantity 



fi{RLj,(f>; Q,)Rduj( 



(32) 



gives the probability that the object, which is at a distance 
R, has angular speed between cj and oj + dco and is moving in 
a direction between 4> and 4>+d<j). This expression is adapted 
for use with our new coordinate system. 

The third factor, C3, is the probability that the object 
triggers another event. This is simply given by 

Combining this with (32) and integrating over V{p) yields 
the quantity C2 x C3: 



C2 X Ca 



LA 
4irR^ 



dw dd> dt. 



(34) 



which illustrates the benefits of our choice of coordinates 
{uj,4>,t)p. In words, C2 x C3 is the probability that an i- 
type object with luminosity L, distance R, and starting at 
the location Clp, moves into the region V{p) and triggers an 
event. 

Now we can put together all the factors which make up 
C(p; V) to find 



00 00 



C{p;V) = pV{p) + 
= pV{p) + 



' L=0 R=0 

A^yi 

AtvJ p 

' L=0 fi=0 V(p) 



Cl X C2 X Cs 



{R,^p) 



X Pi,L{L)L^ fi{Ru,(t>;Q.p)dLod(t>dtdLdR. (35) 

In practice, since resolved objects will be removed from 
the sky map, the lower limit of the R integral should be cut 
off so that these objects are not counted. If the detector can 
resolve any source with flux greater than Fres t hen the lower 
limit on the R integral should be y/L/4nFias- 

Of course, if ni{R,(l) is cut off at a lower limit Rmin 
and Pi^L{L) is cut off at an upper limit Lmax such that 
(imax/47ri?^in) < Frcs uo chauges need to be made to the 
limits of integration in (35) since all i-type objects will be 
unresolved. 



8.5 The form of ^ in 3 dimensions 

Finally, we can substitute (35) into the definition of the 2- 
point function ^ (2) and arrive at an expression for the 2- 
point spacetime correlation function in 3 dimensions, 



^ = 



1 

4-Kp 



EE 

X Pi,L(L) L'^fi{Rw,cl>; Clp)du}dcj>dt 



ni{R, Clp) 
R 



Y,v{p) 



(36) 



V(p) 



Notice that the detector area A has cancelled when using p, 
the average fiux per solid angle, instead of p. It is also ap- 
parent that the contribution to the ^ from different classes 
of objects as well as from objects of difi'crent distances and 
luminosities is additive. The observed 2-point function is 
simply the sum of contributions from different types of ob- 
jects. As expected, the correlation is increased for brighter- 
appearing objects as is seen by the presence of and R~^. 
The interplay between distance and angular speed appears 
in the argument of the velocity distribution 

The expression for ^ given by (36) is a main result of 
this paper. In its general form, however, it is fairly opaque. 
We can get a qualitative feel for the 2-point function by 
calculating ^ for a very simple model whore wo have only 
one class of objects. Those objects have a constant rmmber 
density n and are found only at distances between _Ri and 
R2- The intrinsic lununosity of all the objects will be fixed 
at A so that Pi^L{L) = &{L — A). We choose an isotropic 
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Maxwell-Boltzmann velocity distribution. Projected into 2 
dimensions it becomes the Rayleigh distribution: 



f(v, (j>; dv d(f) - 



dv- 



2tt' 



(37) 



independent of Cl. Finally we choose V{p) to be given by 
(24) in which (j> runs from to 2-k. With these choices no 
quantity in (36) depends on p so both sums over p disappear. 

Let's look at the limiting form for ^ by choosing an 
infinitesimal volume for V{p) where ci;2 = i^i + dij and t2 = 
t\ + dt. Dropping the subscripts on lo\ and t\ we have the 
following expression for the 2-point function, 



C(^,i) 



l-Kt sin(ajt) 
1 

27rt sin(ajt) 

R2IJJ 
V2a 



^ V 

47rp / 



R2 



dRRu 



1 

47r/3 



2 >2 



Erf 



- Erf 



Rioj 
V2a 



(38) 



Note that the contribution to 5(cj,t) from objects at 
different distances serves to smear the influence of f{v) so 
that ^ is not simply proportional to the velocity distribution 
as it was in the 2-dimensional model. There is, however, a 
functional similarity to the 2-dimensinal case: 



^3D 



1 nX 



2 fiR^) 

Pv, (using pi — nX in Eq. 15) 



V{p) p 
1 nX' 

where f{Rio) represents the smeared velocity distribution. 



(39) 
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Figure 7. Results from a simulation of moving objects in the 
solar system, along with stationary sources and random noise. 
The correlation function is plotted for angular velocities between 
and 500° /yr. Red x's represent the theoretical value of ^ calcu- 
lated from (36) while the blue squares show the measured value 
of ^ from the sky map. The width of each angular velocity bin is 
20°/yr. Error bars are derived using (16). The spike at zero angu- 
lar velocity is due to the presence stationary background sources. 
The correlation function is also non-zero between u = 196°/yr 
and 389°/yr, corresponding to moving sources orbiting between 
0.95 and 1.5 AU. 



8.6 Mock Fermi search for solar system bodies 

In order to verify the formulation of ^ in (36) we simulate 
a mock 5-year Fermi observation of nearby moving gamma- 
ray sources. These objects might correspond to a population 
of bodies in the asteroid belt (see Sec. 10 for motivation). 

For simplicity, the detector is a stationary observer at 
the center of the solar system and the moving objects are 
placed on circular orbits with Keplerian velocities deter- 
mined by their distance from the Sun {v cx R~^^^). The ob- 
jects are distributed uniformly in a disk with uniform surface 
density between distances of 0.95 AU and 1.5 AU. For geo- 
metric simplicity, the inclination angles of the orbits are ran- 
dom so that the flux is statistically isotropic. Each object has 
the same luminosity and the closest object (at 0.95 AU) has 
a photon flux of 1.8 x 10~^°cm~^s~^. Note that this flux is 
below the point source detection limit of Fermi so that none 
of these moving objects would be individually identifled as 
localized sources®. The sky contains 7371 objects so that, for 
a 5 year Fermi observation (effective area ~ 2000 cm^), the 
population of moving objects contributes roughly 2.5 x 10^ 
events to the sky map. In addition, as in the 2-dimensional 
simulation, we include a population of stationary objects as 

® In fact, moving sources will be more difficult to detect than sta- 
tionary ones because of their apparent motion, i.e. standard point 
source analysis may be inefficient at detecting moving sources. 



well as completely random events. The stationary sources 
generate detected events at an average rate of 0.2 events per 
year and are distributed isotropically. The stationary and 
random components each comprise about 1.25 x 10^ events 
so that the sky map contains about 5 x 10^ events, 50% 
from moving objects, 25% from stationary sources, and 25% 
random events. 

In computing the correlation function we use spacetime 
volumes V{p) given by (24) with ti — and AT = t2—ti = 
0.015 yr (~ 5.5 days). The angular velocity bins run from 
to 500°/yr in steps of 20°/yr. Results of the measurement 
of ^ for the different angular velocities are shown with blue 
squares in Fig. 7. As with the 2-dimensional simulations, 
the error bars are found using (16). The presence of both 
the stationary and moving sources can be easily seen in the 
shape of the correlation function. 

The theoretical value of ^ based on the properties of the 
sources is a straightforward application of (36). The number 
density of moving objects is 

m = / ^objs/[2^i?(i?i - Rl)] Ri<R<R2 .^Q. 
1 otherwise. 

Here, Ri = 0.95 AU, R2 = 1.5 AU, and Nohj^ = 7371. The 
luminosity function and the velocity distribution are delta 
functions: 

Pl{L) = 5(L-4.45 X lO'^ec"^) (41) 
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fivA) = ^6{v-vo{R/Ror'^^). (42) 

In the above, Ro = 1 AU and vo = 27rJ?o/yr = uJoRo- The 

average event rate p is estimated by dividing the total num- 
ber of events by the sohd angle of the sky map and by the 
observation time: p = 5 x 10"'/(47r x Syr). Carrying through 
the calculation of (36) yields f for the particular choice of 
Vip): 

^/ f LA y 

^(a;i,a;2;AT) = 4^^2y(^^,^^. at) [^J 

^objs ( 1 1\ 
-"-2 ~ ^1 V / 

The quantity V{loi,uj2; AT) is the volume of the spacetime 
region given in (27), Ra = Max{Ri,Ro{uJo/uj2f^^), Rb = 
Min(ii2,-Ro(wo/wi)^/^), and ^ = if Ra > Rb- In (43) the 
quantities Ri, R2, Ra, and Rb are in units of J?o = 1 AU. 

The correlation for the stationary objects is much sim- 
pler. It is equal to zero unless oJi = 0, in which case it is 
given by 

where A^stat = 1.25 x 10®, the number of stationary objects 
and A = 0.2/yr, the detected event rate for each stationary 
object. 

The total correlation function will be the sum of the 
correlation functions for each component. This sum is plot- 
ted in Figure 7 as red x 's, demonstrating that the formalism 
predicts the correct value for the correlation function. 

Although intended as a toy model, this simulation cap- 
tures the essential components of a large area analysis of 
Fermi data. In reality, Fermi has detected far more than 
5 X 10® events. If all the components in our toy model were 
scaled up appropriately the detection of ^ would be even 
more significant. 

8.7 Errors and flux-limited vs. counts-limited 
surveys 

In three dimensions the errors on f given by Eqs. 16 & 17 
also apply. As above, choosing the regions V{p) requires bal- 
ancing a large signal to noise ratio against having many inde- 
pendent choices of V{p). In order to make more independent 
measurements of ^ the size of V{p) must decrease. 

A larger V{p) has its advantages and disadvantages. A 
large volume V will decrease the fluctuations in ^ because 
more events are collected in each such volume (the signal-to- 
noise contains a V^^^ factor). On the other hand, having a 
lot of events in V(j)) which are uncorrelated to the event at 
p will dilute the amplitude of ^ because of the pV{p) term in 
the denominator (c.f. the definition of ^ (2)). There is then 
a tradeoff between the fiuctuations in ^ and the amplitude 
ofe 

If the sky map has a large number of events then it is 
permissible to choose V{p) to be small and still have small 
fiuctuations in ^. In the opposite limit, if the sky map is 
"counts-limited" then it will be necessary to choose V{p) to 
have a large volume. The safest method for deciding is to 
run realistic simulations for various combinations of physical 
parameters and experiment with different choices for V{p). 



There is an additional requirement on V{p) which de- 
pends on the detector's resolution. If one chooses V(j)) to be 
very small (in the angular sense) then one is essentially ask- 
ing the detector to distinguish events at this angular scale. 
The detector has a smallest "pixel size" and V{p) cannot be 
smaller than that. 

The most convenient choice for V{p) when calculating ^ 
according to (36) is given by (24). Unfortunately, this choice 
is inconvenient when dealing with a detector with a finite 
angular resolution (a real detector). The projection of the 
spacetime region V{uji,uj2;ti, ^2) onto the celestial sphere 
must have an angular size no smaller than the detector's 
angular resolution. However, for fixed Aw = 0)2 — (jJi and 
At = t2 — ti, changing cui and ti will change the projected 
angular size of V{p). The bin sizes Alo and At must be 
varied with wi and ti . An estimate of this constraint is that 
the angular resolution of the detector be no worse than 9 « 
A(Lut) ~ CuAt + Aujt, where Auj = iU2 — 1^1 and At — t2 —t2. 

All of these choices are part of the auah'sis, not the 
collection, of the data. If the diffuse background events are 
already in hand one can experiment with different choices 
for the V{pys to find the right balance between signal-to- 
noise and number of independent measurements of ^ while 
maintaining the detector resolution constraint. 

Of course, it is best to use PSF information instead of 
an assumption of "pixel size" . We have discussed this option 
for the 2-dimensional case (see (18)). The generalization to 
3 dimensions is straightforward. 



9 GENERALIZATIONS 

There are several ways to make this technique more power- 
ful. Here we mention two: the inclusion of spectral data and 
the use of n-point functions. 

9.1 Including spectral information 

Not only do sky surveys keep track the direction and time of 

each photon they receive, they can also measure wavelength 
(or energy of the photon) . The easiest way to make use of 
this information is to note that the above analysis holds 
for every wavelength separately. One can bin the events by 
energy, make separate sky maps for each energy bin, and 
then compute the 2-point function for each of the maps. 
Typically, this procedure will add more data points than 
free parameters: the same distributions Ui and fi are used 
for different energy bins. Only the luminosity functions will 
vary, though the physical parameters in Pi^L are likely to 
be universal over all energy bins. Thus, ^ measured at one 
energy will be related to ^ measured at another. As a re- 
sult, an analysis which includes event energies can help in 
untangling the different components of the background. 

9.2 n-point functions 

Another generalization of the 2-point function is, naturally, 
the n-point function. One asks, "Given an event at p what 
is the probability of finding events in Vi{p) and in V2(p)?" 
If the objects move in straight lines then this probability 
will spike when p, Vi{p), and V2(p) lie along a straight line. 
The jump from n = 2 to n = 3 is significant for this reason 
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— every pair of points is coUincar but not every trio. The 
downside of measuring n-point functions (besides the com- 
putational cost) is that they require a much larger number 
of events to overcome statistical fluctuations. Recall that in 
our 2-dimensional toy model only 5% of events came from 
objects which generated more than one event and that of 
these events, 95% came from objects which generated ex- 
actly two events. Therefore, only 0.25% of the events in 
the map came from objects which generated three or more 
events. Although they will be slightly more cumbersome, 
analytic forms for these higher correlation functions can be 
found by applying the same reasoning we used for the 2- 
point function. 



10 DISCUSSION AND CONCLUSIONS 

We present a new tool, based on the familiar 2-point corre- 
lation function, which can be applied to astrophysical maps 
of diffuse emission. The measured quantity ^ is designed to 
detect the presence of moving objects, each of which is too 
dim to be resolved individually. We derived the form of ^ 
based on the physical parameters which describe the classes 
of objects which might be present in the sky (36). A mea- 
surement of ^ along with the theoretical prediction for ^ can 
be used to find best-fit quantities for the physical parameters 
describing the populations of objects. We emphasize that all 
the technology invented to study the angular 2-point corre- 
lation function can be directly applied to the generalization 
to the spacetime 2-point correlation function. 

There are numerous applications of the derived formal- 
ism. An obvious place to start is the diffuse gamma-ray back- 
ground measured by the Fermi-LAT instrument. The all-sky 
capabilities of LAT, coupled with its high angular resolu- 
tion provide a convenient testbed where this technique can 
be applied. The interesting question is what kind of sources 
contribute to the gamma-ray background and also exhibit 
proper motion over the duration of observation. 

One potential source is the generation of gamma-rays 
from cosmic-ray interactions in rocky debris present in the 
solar system. Cosmic ray interactions with nuclei on a solar 
system body lead to hadronization, and the subsequent de- 
cay of neutral pions to a photon final state (Stecker 1970; 
Stephens & Badhwar 1981; Dermer 1986; Kelner, Aharonian 
& Bugayov 2006). A detection of a large population of these 
sources is important as it provides information about the 
origin of the solar system and its evolution with time, as 
well as the energy spectrum and composition of the incident 
cosmic ray flux. 

The detection of gamma-rays from cosmic ray interac- 
tions with solar system bodies has been discussed in the 
context of past measurements by the Energetic Gamma Ray 
Experiment Telescope (EGRET) on board the Compton 
Gamma-ray Observatory, and measurements with Fermi- 
LAT (Moskalenko et al. 2008; Moskalenko & Porter 2009; 
Giglietto et al. 2009c). Sources include small objects in the 
main asteroid belt, Trans- Neptunian objects in the Kuiper 
belt, as well as objects in the Oort cloud, including icy bod- 
ies such as comets. It was shown that for objects where the 
cosmic ray cascade fully develops (objects with size greater 
than ~ 1 m) it may be possible for Fermi to detect the 
cumulative gamma-ray emission from a collection of such 



bodies. These estimates are based on the distribution and 
composition of objects. Even though both of these quanti- 
ties are partially constrained for objects in the main asteroid 
belt, large uncertainties arc present for the populations in 
the Kuiper belt and the even more speculative Oort cloud. 
It is conceivable that a very large number of bodies may be 
present in the outskirts of the solar system. 

The proximity of these populations makes them ideal 
for an application of the spacetime correlation function, as 
each source will traverse an angular distance which is larger 
than the angular resolution limit of Fermi. Typical angular 
displacement (assuming Keplerian orbits) of an object at 
distance d from the Sun is 9 = 27r rad(Ar/yr)(d/AU)"^/^ 
during the course of an integration for time AT. The com- 
position of these objects can be assumed to be similar to 
the composition of the Moon, though their mass density 
varies considerably. This similarity in composition is con- 
venient as the gamma-ray flux due to cosmic interactions 
with the lunar rock is well understood (Moskalenko & Porter 
2007; Morris 1984) (see also (Giglietto et al. 2009a,b)). If 
we assume that the spectral shape of the gamma-ray emis- 
sion from solar system bodies is similar to that of the rim 
of the Moon (emission above 600 MeV is dominated by 
the rim of the Moon rather than the lunar disc) and we 
scale the flux from the object to the flux from the Moon 
($M = 1.1 X 10"®cm"^s"\ Giglietto et al. (2009a)), the 
flux from an object of radius r at distance d would then 
be $ = ^M{r/rM)(dM /dY . For a distance to the Moon of 
dm = 0.0024 AU and a lunar radius of vm = 1740 km, the 
total number of photons per year detected by the Fermi- 
LAT instrument (with an orbit-averaged effective area of 
2000 cm^) is $ ^ 2 X 10"''yr-i(r/km)(rf/lAU)"^ There- 
fore, given this information, one can apply the spacetime 
correlation function to determine the abundance and radial 
distribution of solar system objects that contribute to the 
gamma-ray background (Geringcr-Sameth & Koushiappas 
2011). It is important to note that even though a theoret- 
ical estimate of ^ requires knowledge of the objects one is 
searching for, the measurement of ^ requires no such knowl- 
edge. 

Similar arguments can be used in search of the ener- 
getic neutrino signal from cosmic ray interactions with solar 
system bodies. The decay of kaons to charged pions leads to 
an energetic signal with a spectral signature that is different 
from the cosmic ray neutrino flux expected from spallation 
of nuclei. Therefore, energetic neutrinos from cosmic ray in- 
teractions with solar system bodies should be present in the 
signal measured by IceCube (Icecube Collaboration et al. 
2006). The sources of these neutrinos will traverse an an- 
gular distance based on the distance of the source from the 
Sun, and therefore the spacetime correlation function de- 
rived here can be used in search of these sources. However, 
as in the case of gamma-rays, the uncertainties in the distri- 
bution and composition of small solar system bodies make 
predictions for such signal difficult. Nevertheless, a blind 
analysis of neutrino events from IceCube could place con- 
straints on the parameters that describe the different popu- 
lations of small bodies in the solar system. 

Another application is in the search for primordial black 
holes in the solar neighborhood. Primordial black holes may 
form in the early Universe through the collapse of large pri- 
mordial fluctuations (Hawking 1971). Current bounds on the 
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abundance of such black holes are of order S^pbh ~ 10^ for 
most of the range of black hole masses (Lacki & Beacom 
2010). If primordial black holes exist in an otherwise dark 
matter dominated Universe, they will acquire a dark matter 
halo (Mack, Ostriker & Ricotti 2007; Ricotti 2007). Dark 
matter annihilation around primordial black holes and/or 
high density ultracompact halos will result in gammarray 
emission (Ricotti & Gould 2009; Scott & Sivertsson 2009). 
Such objects with very small mass will in fact be very dense 
and survive in the Milky Way halo. If we assume that pri- 
mordial black holes trace the distribution of dark matter in 
the Milky Way we can use their abundance to determine the 
angular distance that a black hole may traverse in a given 
time interval. For simplicity, let's assume that primordial 
black holes have mass MpBH = lO-i^Mo, QpBH = 10-^ 
and that the local dark matter density is 0.01 Mqpc~^. 
Then the mean distance between primordial black holes in 
the solar neighborhood is ~ 10~^pc. Assuming that this is 
the maximum distance to a primordial black hole, and that 
the mean velocity of primordial black holes is similar to the 
mean velocity of dark matter, i.e., 220 km/s, then the an- 
gular displacement of these gamma ray sources can be as 
large as 4.5 degrees in 10 years. As the angular resolution of 
Fermi is significantly less for energies greater than 1 GeV, 
constraints on the abundance and size of these black holes 
can be placed by applying the spacetime correlation function 
to the LAT all-sky map. 

A more speculative contribution to the gamma-ray 
background is from dark matter halos formed on scales close 
to the cutoff scale of the dark matter power spectrum. These 
objects typically have sub-solar masses (Schmid, Schwarz & 
Widerin 1999; Hofmann, Schwarz & Stocker 2001; Green, 
Hofmann & Schwarz 2004, 2005; Loeb & Zaldarriaga 2005; 
Chen, Kamionkowski & Zhang 2001; Profumo, Sigurdson & 
Kamionkowski 2006). Even though their survival and abun- 
dance in the present-day Milky Way halo is unknown, it is 
possible that dark matter annihilation in these high-density 
objects may contribute to the gamma-ray background (Fieri, 
Bertone & Branchini 2008; Ando et al. 2008). The probabil- 
ity that such sources will exhibit spatial motion in the du- 
ration of the Fermi-LAT mission is directly linked to their 
abundance, and thus the use of the correlation function can 
provide information on the survival rate of these extremely 
early-forming objects. 

The spacetime correlation function can be applied to 
lensing surveys to search for compact objects in the Milky 
Way. Past studies suggest that up to 20% of unseen matter 
is in the form of Massive Compact Halo Objects (MACHOs) 
(Alcock et al. 2000; Uglcsich et al. 2004). With the advent 
of dedicated surveys e.g., LSST, (LSST Science Collabora- 
tions et al. 2009), as well as astrometric missions such as 
SIM (Unwin et al. 2008) and Gaia (Lindegren et al. 2008), 
it will be possible to generate time-domain maps of lens- 
ing events in dense stellar fields. Such information can be 
used to probe correlated events originating from the spatial 
translation of compact objects, thus probing the projected 
velocity distribution of the compact population in the Milky 
Way. In addition, it may also be possible to place constraints 
on the density, abundance and distribution of dark matter 
substructure (Erickcek &: Law 2011). 

Throughout the development of the analysis we as- 
sumed that the event rate due to any source was constant in 



time. There are many classes of astrophysical objects with 
time-dependent emission. Most notably, unresolved pulsars 
are thought to contribute to the diffuse gamma-ray back- 
ground (e.g. Watters & Romani (2011); Fauchcr-Gigucrc & 
Loeb (2010)). While these sources will not exhibit proper 
motion over the course of observations the temporal corre- 
lations of their emitted photons may be discovered through 
techniques based on the ones presented here (Geringer- 
Sameth & Koushiappas 2012). Essentially, one chooses the 
volumes V{p) according to (3) (illustrated in the left panel 
of Fig. 3), but with a non-trivial slicing along the time axis. 
Such a V{p) picks up on stationary objects which exhibit 
correlations within their photon time series. 

The power of this analysis for untangling the contri- 
bution of different classes of sources requires that each class 
have "different enough" velocity, luminosity, and spatial dis- 
tributions. For example, if two classes have similar velocity 
and spatial distributions then one may as well just treat 
them as a single class with a modified luminosity function. 
This points to a problem that is likely to be encountered 
in many realistic astrophysical applications: the angular ve- 
locities of almost all objects will be much too small to be 
resolved by a detector. That is, when one combines the ve- 
locity distribution /» with the spatial distribution m in (36) 
it may be that ^ = at all angular velocities except in a tiny 
range near u = 0. This is because virtually all of the objects 
have distances and speeds such that their apparent proper 
motion is below the angular resolution of the detector. A 
large degeneracy is created and it will be impossible to pull 
out information about any specific class of objects. The fact 
that ^ is not zero at a; = indicates the existence of ob- 
jects. However, without being able to measure the shape of 
^ for different angular speeds cj the 2-point function loses its 
value as a tool to untangle the contributions from different 
classes of objects. 

Of course, as the resolutions of detectors improve, the 
2-point function becomes more useful. It is a straightforward 
task to calculate ^(wi , 0)2 ; i i , ^2 ) for specific classes of objects 
and find out over what ranges of u and t the correlation 
drops to zero. For example, if ^ goes to zero around (wi = 
ui\ ti = t') then a detector which has resolution better than 
6 ~ 10' t' can measure the shape of ^{iO,t) as it goes from a 
maximum at (wi = 0,ti — 0) to zero at {uji = Lo',ti — t'). 

In summary, we introduced the spacetime correlation 
function, a statistical tool that can be used to search for 
the presence of moving, flux-unresolved sources in a dif- 
fuse background. This formalism has numerous applications. 
With large area sky surveys and long duration baselines the 
spacetime correlation function can be used to disentangle 
the contributions from spatially moving sources, and may 
aid in the discovery of new sources. 

We thank the anonymous referee for comments and 
suggestions that improved the quality of this manuscript. 
We acknowledge useful conversations with John Beacom, 

.Jacqueline Chen, Richard Cook, Ian DcH'Antonio, Scott Do- 
delson, Andrew Favaloro, Salman Habib, Arthur Kosowsky, 
David Laidlaw, Miguel Morales, Louie Strigari, Andrew 
Zentner. SMK and AGS are funded by NSF PHYS-0969853 
and by Brown University. 



© 0000 RAS, MNRAS 000, 000-000 



Detecting unresolved moving sources in a diffuse background 17 



REFERENCES 

Abdo A. A. et al., 2009, Physical Review Letters, 103, 
251101 

Alcock C, Allsman R. A., Alves D. R., Axelrod T. S., 

Becker A. C, Bennett D. P., Cook K. H., et al., 2000, 

Astrophys. J., 542, 281 
Ando S., 2009, Phys. Rev., D80, 023520 
Ando S., Kamionkowski M., Lee S. K., Koushiappas S. M., 

2008, Phys. Rev. D, 78, 101301 
Ando S., Komatsu E., Narumoto T., Totani T., 2007, Phys. 

Rev., D75, 063519 
Ando S., Pavlidou V., 2009, Mon. Not. R. Astron. Soc, 

400, 2122 

Barrow J. D., Bhavsar S. P., Sonoda D. H., 1984, Mon. 
Not. R. Astron. Soc, 210, 19P 

Baxter E. J., Dodelson S., Koushiappas S. M., Strigari 
L. E., 2010, Phys. Rev. D, 82, 123511 

Bernardi G., Carretti E., Fabbri R., Sbarra C, Cortiglioni 
S., 2005, Mon.Not.RoyAstron.Soc.Lett., 364, L5 

Bernstein G. M., 1994, Astrophys. J., 424, 569 

Bhattachaxya D., Sreekumar P., Mukherjee R., 2009, Re- 
search in Astronomy and Astrophysics, 9, 1205 

Chen X., Kamionkowski M., Zhang X., 2001, Phys. Rev. 
D, 64, 021302 

Chiang J., Mukherjee R., 1998, Astrophys. J., 496, 752 
Colafrancesco S., Blasi P., 1998, Astroparticle Physics, 9, 
227 

Cruz M., Viclva P., Martmcz-Gonzalcz E., Barreiro R. B., 
2011, Mon. Not. R. Astron. Soc, 412, 2383 
Cuoco A., Brandbyge J., Hannestad S., Haugboelle T., 

Miele G., 2008, Phys. Rev., D77, 123518 
Cuoco A., Hannestad S., Haugb0llc T., Miclc G., Scrpico 
P. D., Tu H., 2007, J. Cosmology Astropart. Phys., 4, 13 
Dermer C. D., 1986, Astron. & Astrophys., 157, 223 
Dobler G., Draine B. T., Finkbeiner D. P., 2009, Astro- 
phys.,!., 699, 1374 
Dobler G., Finkbeiner D. P., 2008, Astrophys. J., 680, 1222 
Dodelson S., Belikov A. V., Hooper D., Serpico P., 2009, 

Phys.Rev., D80, 083504 
Erickcek A. L., Law N. M., 2011, Astrophys. J., 729, 49 
Erikson H. K., Lilje P. B., Banday A. J., Gorski K. M., 

2004, ApJS, 151, 1 
Fall S. M., Tremaine S., 1977, Astrophys. J., 216, 682 
Fauchcr-Giguere C.-A., Loeb A., 2010, J. Cosmology As- 
tropart. Phys., 1, 5 
Fisher K. B., Davis M., Strauss M. A., Yahil A., Huchra 

J., 1994, Mon. Not. R. Astron. Soc, 266, 50 
Fornasa M., Pieri L., Bertone G., Branchini E., 2009, Phys. 

Rev., D80, 023518 
Fry J. N., Gaztanaga E., 1993, Astrophys. J., 413, 447 
Geringer-Sameth A., Koushiappas S. M., 2011, in prepara- 
tion 

— , 2012, Mon. Not. R. Astron. Soc, 421, 1813 
Giglictto N., et al., 2009a, ArXiv c-prints 0907.0543 
— , 2009b, ArXiv e-prints 0912.3734 
— , 2009c, ArXiv e-prints 0907.0541 
Gold B. ct al., 2009, ApJS, 180, 265 

Green A. M., Hofmann S., Schwarz D. J., 2004, Mon. Not. 

R. Astron. Soc, 353, L23 
— , 2005, J. Cosmology Astropart. Phys., 8, 3 
Hamilton A. J. S., 1993, Astrophys. J., 417, 19 



Hawking S., 1971, Mon. Not. R. Astron. Soc, 152, 75 
Hofmann S., Schwarz D. J., Stocker H., 2001, Phys. Rev. 

D, 64, 083507 
Hooper D., Scrpico P. D., 2007, JCAP, 0706, 013 
Icecube Collaboration, Achterberg A. et al., 2006, As- 
troparticle Physics, 26, 155 
Inoue Y., Totani T., 2009, Astrophys. J., 702, 523 
Kelner S. R., Aharonian F. A., Bugayov V. V., 2006, Phys. 

Rev. D, 74, 034018 
Komatsu E. et al., 2009, ApJS, 180, 330 
Koushiappas S. M., 2006, Physical Review Letters, 97, 
191301 

— , 2009, New Journal of Physics, 11, 105012 
Lacki B. C, Beacom J. F., 2010, Astrophys. J. Lett., 720, 
L67 

Landy S. D., Szalay A. S., 1993, Astrophys. J., 412, 64 
Lee S. K., Ando S., Kamionkowski M., 2009, J. Cosmology 

Astropart. Phys., 7, 7 
Lindegren L. et al., 2008, in lAU Symposium, Vol. 248, lAU 

Symposium, W. J. Jin, I. Platais, &: M. A. C. Ferryman, 

ed., pp. 217-223 
Ling E. N., Barrow J. D., Frenk C. S., 1986, Mon. Not. R. 

Astron. Soc, 223, 21P 
Loeb A., Waxman E., 2000, Nature, 405, 156 
Loeb A., Zaldarriaga M., 2005, Phys. Rev. D, 71, 103520 
LSST Science Collaborations et al., 2009, ArXiv e-prints 

0912.0201 

Lupton R., 1993, Statistics in theory and practice, Lupton, 
R., ed. 

Mack K. J., Ostriker J. P., Ricotti M., 2007, Astrophys.J., 

665, 1277 

Miniati F., 2002, Mon. Not. R. Astron. Soc, 337, 199 
Miniati F., Koushiappas S. M., Di Matteo T., 2007, Astro- 
phys. J. Lett., 667, LI 
Mo H. J., Jing Y. P., Boerner G., 1992, Astrophys. J., 392, 
452 

Moore A. W. et al., 2001, in Mining the Sky, A. J. Banday, 
S. Zaroubi, & M. Bartelmann, ed., p. 71 

Morales M. F., Williams D. A., De Young T., 2004, As- 
troparticle Physics, 20, 485 

Morris D. J., 1984, J. Geophys. Phys., 89, 10685 

Moskalenko I. V., Porter T. A., 2007, Astrophys. J., 670, 
1467 

— , 2009, Astrophys. J. Lett., 692, L54 

Moskalenko I. V., Porter T. A., Digol S. W., Michelson 

P. F., Ormes J. F., 2008, Astrophys. J., 681, 1708 
Miicke A., Pohl M., 2000, Mon. Not. R. Astron. Soc, 312, 

177 

Narumoto T., Totani T., 2006, Astrophys. J., 643, 81 
Padovani P., Ghisclhni G., Fabian A. C, Celotti A., 1993, 

Mon. Not. R. Astron. Soc, 260, L21 
Pavlidou v.. Fields B. D., 2002, Astrophys. J. Lett., 575, 

L5 

Peebles P. J. E., 1973, Astrophys. J., 185, 413 

Pieri L., Bertone G., Branchini E., 2008, Mon. Not. R. 

Astron. Soc, 384, 1627 
Profumo S., Sigurdson K., Kamionkowski M., 2006, Phys- 
ical Review Letters, 97, 031301 
Ricotti M., 2007, Astrophys. J., 662, 53 
Ricotti M., Gould A., 2009, Astrophys. J., 707, 979 
Schmid C, Schwarz D. J., Widerin P., 1999, Phys. Rev. D, 
59, 043517 



© 0000 RAS, MNRAS 000, 000-000 



18 Alex Geringer-Sameth and Savvas M. Koushiappas 



Scott p., Sivertsson S., 2009, Physical Review Letters, 103, 
211301 

Siegal-Gaskins J. M., 2008, JCAP, 0810, 040 

Sicgal-Gaskins J. M., 2010, in Astronomical Society of the 
Pacific Conference Scries, Vol. 426, 2009 Snowbird Parti- 
cle Astrophysics and Cosmology Workshop (SNOWPAC 
2009), D. B. Kieda & P. Gondolo, ed., p. 87 

Siegal-Gaskins J. M., Pavlidou V., 2009, Pliys. Rev. Lett., 
102, 241301 

Stecker F. W., 1970, Astr. and Space Sci., 6, 377 
Stecker F. W., Salamon M. H., 1996, Astrophys. J., 464, 
600 

Stephens S. A., Badhwar G. D., 1981, Astr. and Space Sci., 
76, 213 

Szapudi I., Colombi S., 1996, Astrophys. J., 470, 131 
Szapudi S., Szalay A. S., 1998, Astrophys. J. Lett., 494, 

L41+ 

Taoso M., Ando S., Bcrtone G., Profumo S., 2009, Pliys. 

Rev., D79, 043521 
Uglesich R. R., Crotts A. P. S., Baltz E. A., de Jong J., 

Boyle R. P., Corbally C. J., 2004, Astrophys. J., 612, 877 
Unwin S. C. et al., 2008, Publ. Astron. Soc. Pac., 120, 38 
Venters T. M., 2010, Astrophys. J., 710, 1530 
Watters K. P., Romani R. W., 2011, Astrophys. J., 727, 

123 

Zehavi I. et al., 2002, Astrophys. J., 571, 172 

Zhang L. L., Pen U.-L., 2005, New Astronomy, 10, 569 



© 0000 RAS, MNRAS 000, 000-000 



